Critical dynamics of the simple-cubic Heisenberg antiferromagnet RbMnF 3 : 

Extrapolation to q=0 
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Monte Carlo and spin dynamics simulations have been used to study the dynamic critical behav- 
ior of RbMnF3, treated as a classical Heisenberg antiferromagnet on a simple cubic lattice. In an 
attempt to understand the difference in the value of the dynamic critical exponent z between exper- 
iment and theory, we have used larger lattice sizes than in our previous simulations to better probe 
the asymptotic critical region in momentum. We estimate z — 1.49 ± 0.03, in good agreement with 
the renormalization-group theory and dynamic scaling predictions. In addition, the central peak 
in the dynamic structure factor at T c , seen in experiments and previous simulations, but absent 
in the renomalization-group and mode-coupling theories, is shown to be solely in the longitudinal 
component. 



I. INTRODUCTION 

The dynamic critical behavior of RbMnF3 , a close real- 
ization of the isotropic, simple cubic Heisenberg antifer- 
romagnet, has been the subject of several experimental 
(see Coldea et ali and references therein), and theoretical 
studiesS*^*^. The real time dynamics of this system is 
governed by coupled equations of motion for the magnetic 
ions. In the classification of Hohenberg and Halperin 5 , 
the critical dynamics of this system pertains to class G, 
for which the order parameter (the staggered magnetiza- 
tion) is not conserved. 

Dynamic critical behavior can be characterized by a 
dynamic critical exponent z, and for RbMnFa the most 
precise experimental estimate— is z = 1.43 ± 0.04. This 
estimate is slightly below the predicted value^^ of z = 
1.5 for an isotropic three-dimensional Heisenberg antifer- 
romagnet. Our previous estimate using spin-dynamics 
simulation^ (obtained before we knew the latest experi- 
mental result) was z = 1.43 ± 0.03 in good agreement 
with the experimental value. The slight disagreement 
between theory and both the experiment and the simu- 
lation was perplexing, but one possible explanation was 
that neither of these latter studies probed the asymp- 
totic critical region sufficiently far. To check this, neu- 
tron scattering experiments would have to be performed 
with smaller momentum transfer q, a task that can be 
quite challenging. We can also resort to spin dynamics 
simulations of the model on larger lattice sizes which al- 
low access to smaller q values. 

Spin dynamics simulations 8 have proven to be an ef- 
fective tool to study dynamic behavior of magnetic sys- 
tems. Direct and quantitative comparisons of magnetic 
excitation dispersion curves and dynamic structure fac- 
tor line shapes from experiments! and spin dynamics 
simulations have shown good agreement, with no ad- 
justable parameters^. At the critical temperature T c , 
both experiment and simulation find that the average dy- 
namic structure factor S(q, u>), for momentum q and fre- 
quency w, has a spin- wave peak and an additional central 



peak not predicted by either the renormalization-group 
theory^ or mode-coupling theory^. In a polarized neu- 
tron scattering experiments below T c , the quasi-elastic 
peak has been shown to be longitudinal in character. It 
is thus interesting to investigate whether the central peak 
at T c originates in the longitudinal, in the transverse or in 
both components with respect to the staggered magneti- 
zation. Another motivation for separating the longitudi- 
nal and transverse components of the dynamic structure 
factor is to test the theoretical prediction^ that below T c 
the longitudinal momentum-dependent susceptibility be- 
haves as x L (q) ~ whereas the g-dependence of the 
transverse component is x T (<?) ~ V? 2 - Although exper- 
imental measurements^ are consistent with x L (l) ~ 
the lack of reliable small wave vector measurements hin- 
dered a conclusive experimental test of this predicted 
divergence. The experimental value for the transverse 
component is x T (g) ~ i/g 1 - 91 ^ - 05 . 

In this paper we use spin dynamics simulations to 
study the dynamic critical behavior of the isotropic 
Heisenberg antiferromagnet on the simple cubic lattice, 
using larger lattices than in our previous simulations, as 
motivated above. We investigate the longitudinal and 
transverse components of the dynamic structure factor, 
and compare the g-dependence of the susceptibility with 
the predicted forms both at and below T c . For brevity, 
we do not present our methodology in great detail here; 
it has been described in earlier work££. 



II. MODEL AND METHODS 

We consider three-dimensional classical spins S r of unit 
length, defined on the sites r of L x L x L simple cubic 
lattices. The interaction is described by a model Hamil- 
tonian written as 



n — j ' S r • S r 



(1) 



where the summation is over pairs of nearest-neighbor 
spins, and J > is an antiferromagnetic exchange 
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coupling. An earlier high-resolution Monte Carlo 
simulationii determined T c = 1.442929(77) J, so any un- 
certainty in the location of the critical temperature of 
this model is negligible. 

The dynamics of the spins is governed by coupled 
equations of motion^ and the dynamic structure factor 
£ (q, lo) is given by the Fourier transform of the space- 
and time-displaced correlation functions C k (r — r',t) = 
(S k (t)S k (0)) - (S*(t))<S*(a)), where k represents the 
spin components. The coupled equations of motion 
were integrated using an algorithm 12 based on 4th-order 
Suzuki- Trotter decompositions of exponential operators, 
with a time step dt. The integrations were performed up 
to time t max , starting from equilibrium configurations 
at temperature T generated with a hybrid Monte Carlo 
methodi. The correlations C k (r — r',t) were computed 
for time displacements ranging from to t cu tof f and the 
canonical ensemble was established by averaging results 
from N different initial equilibrium configurations. We 
used periodic boundary conditions in space and thus we 
could only access a set of discrete values of momentum 
transfer given by q — 2im q /L, where n q = 1,2, ...,L/2. 
Because of computer memory limitations we restricted 
our data to q = (g, 0,0), (q, q, 0), and (q,q,q), which 
correspond to the [100], [110], and [111] directions, re- 
spectively. 

In an attempt to probe the true asymptotic critical 
region in momentum and thus provide a more accurate 
estimate of the dynamic critical exponent, we extended 
our Monte Carlo^ and spin-dynamics simulations^ to 
systems as large as L = 72 at T c . For L — 72, each 
equilibrium configuration was generated with 2500 hy- 
brid Monte Carlo steps, each of which consists of two 
Metropolis sweeps through the lattice and eight overre- 
laxation steps. We used dt = 0.2/J, t max = 1080/ J, 
tcutof / = 1000/ J, and N = 1000. We have also improved 
the statistics of our previous simulations for L = 48 and 
60 by increasing the number of initial configurations to 
N = 1000. With L = 72, the smallest wave vector that 
we can access is q = 2tt/72 = 27r(0.0139), whereas, in 
our units, the smallest q value probed by experiment 1 
was q = 2tt(0.02). 

The dynamic critical exponent z can be ex- 
tracted using dynamic finite-size scaling^, which yields 
^ fc (q^)/x fc (q) = G(ioL*,qL,5 w L*) and O* = 
L z VL k (qL,8 u) L z ), where k represents the polarization, 
S k (q, lo) is the dynamic structure factor convoluted with 
a Gaussian resolution function with parameter 8 U , x fe (q) 
is the total integrated intensity, and is a characteristic 

frequency, defined by f_£ k S k (q,Lo)dLo/2ir — x fc (q)/2. If 
we set 8u = the exponent z can be obtained from the 
slope of a graph of ln(w^) vs ln(L), at fixed qL. To 
test the robustness of the estimate, we have also used 
<L = 0.005(72/L) z , and determined z iterativelji&. As 
in our previous workiS, the dynamic critical exponent is 
estimated using the average dynamic structure factor. 

After understanding the difference in the dynamic crit- 
ical exponent between experiment and theory, we can ex- 



amine other unresolved issues, such as the nature of the 
central peak at T c observed in experiments, but absent 
in theoretical studies. To this end, we analyzed longi- 
tudinal and transverse motions of the spins with respect 
to the staggered magnetization. The dynamics of the 
isotropic Heisenberg model, given by the coupled equa- 
tions of motion, conserves both the total energy and the 
uniform magnetization, which is the order parameter for 
the Heisenberg ferromagnet. However, for the antiferro- 
magnet the order parameter (staggered magnetization) 
is not a constant of the motion; therefore, separating 
components of the spin parallel (longitudinal component) 
and perpendicular (transverse component) to the order 
parameter is challenging. Our approach to determine 
the individual components of the spin motion, and thus 
of 5(q, lo), was to rotate the frame of reference to align 
the z-axis parallel to the staggered magnetization before 
starting the time integrations, to make the z-axis coin- 
cide with the longitudinal direction. As we integrated 
the equations of motion, the direction of the staggered 
magnetization changed slightly because it is not a con- 
served quantity. Therefore, after each integration step 
we re-rotated the frame of reference to re-align the z- 
axis with the staggered magnetization, thereby restoring 
the z-axis as the longitudinal direction. In this part of 
our simulations we used L = 24 with t max — 480/ J, 
t cutof f = 400/ J, dt = 0.2/J, and N = 12,000, in ad- 
dition to L = 36, 48, and 60, with t max = 880/ J, 
t cutof f = 800/ J, dt = 0.2/J, and N = 11,000, 7,000, 
and 3, 000, respectively. We denote the longitudinal 
and transverse components of S(q,u>) as 5 i (q, ui) and 
5 T (q, lo), respectively. 

We have also investigated x L (l) an d X T ( < ?)i the inte- 
grated intensities of S L (q, ui) and S T (q, w), respectively. 
When the decay of the dynamic structure factor is slow, 
as is the case of S T (q, ui) at T c , computing the inte- 
grated intensity is complicated by the fact that at high 
frequencies \^/{2dt) < lo < -k /dt] the approximation of 
the integral Fourier transform as a discrete sum is highly 
dependent on the summation method (e.g. direct sum, 
trapezoidal rule, Simpson's rule, etc.). In terms of cpu 
time usage, it is more efficient to integrate the equations 
of motion with the largest time step dt that still yields a 
stable method and accurate time evolutions; however, the 
maximum accessible frequency is inversely proportional 
to dt. Simulations for L — 24 at T c using a smaller time 
step (dt = 0.1/ J) showed that both the Simpson's rule 
and the trapezoidal rule produce consistent results for 
5 L (q, lo) and S T (q, lo) up to lo ks 15 J. A direct compar- 
ison between the dynamic structure factor obtained with 
dt = 0.1/ J and that obtained with the Simpson's rule 
and the trapezoidal rule using dt = 0.2/J up to lo « 15 J 
shows that the trapezoidal rule gives a better approxima- 
tion for S L (q,co) and S* T (q, lo) at high frequencies. Esti- 
mates of x L (q) an d X T (l) are not significantly dependent 
on the methods used to integrate S L (q, lo) and i5 T (q, lo), 
respectively, and we use Simpson's rule for these integra- 
tions. 
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Our approach to determine x T (q) & t was to compute 
5 T (q, uj) with the trapezoidal rule and then to integrate 
it from ui — to 8 J. Although the intensity of 5 T (q, w) 
has not decayed to zero at uj — 8 J, it is a small fraction 
of the intensity of the spin-wave peak. Data for uj = 4J 
to 8J were fitted with an exponential function, which 
was then used in the integration to w = oo. (An expo- 
nential function was chosen for this fitting because it is 
the simplest one that yielded a good fitting in the range 
of frequency used.) To estimate the longitudinal compo- 
nent X i (<z) at T c we simply used the trapezoidal rule to 
determine S L (q,ui) and then integrated it from uj = 
to lu = 10J, without further high-frequency corrections, 
because 5 L (q, uj) decays quickly in frequency. 

The behavior of x T (l) an d X L (l) below T c is studied 
with spin dynamics simulations at T = 0.5T C using dt — 
0.2/ J and L = 24 and 36, with N = 1000 for each lattice 
size. At this temperature both S L (q, to) and 5 T (q, uj) 
decay quickly, hence these line shapes were obtained with 
the trapezoidal rule and they were then integrated from 
w = 0tow = 10J. 

For comparison, we have also computed the longitudi- 
nal and transverse components of <S(q, uj) with respect to 
the residual uniform magnetization. In this case, the ro- 
tation of the frame of reference to align one axis parallel 
to the uniform magnetization was only required once, be- 
fore starting the time integrations, because the uniform 
magnetization vector is conserved. 



III. RESULTS 

Our results for the average dynamic structure factor 
5(q, uj) show a spin wave and a central peak for T <T C . 
The spin wave dispersion curve for L = 72, T = T c and 
q vectors in the [100] direction, shown in Fig^ was ob- 
tained by fitting 5(q, uj) with Lorentzian spin- wave cre- 
ation and annihilation peaks centered at u) — ±uj s , and a 
Lorentzian peak at uj = 0. We have obtained reasonable 
fittings for q values corresponding to n q = 1 to 7; for 
larger n q , the spin- wave line shapes were not Lorentzian 
and the positions of the spin- wave peaks uj s were read off 
directly without any fitting. We remind the reader that 
the first Brillouin zone edge in the [100] direction occurs 
at n q = L/2. An estimate of the dynamic critical expo- 
nent z is obtained by fitting the dispersion curve with a 
function^ u> s = Dn z q . Since the asymptotic critical region 
corresponds to small values of g, we have used n q — 1 
to 9 in the fitting (dashed line in Fig^l, and obtained 
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z = 1.35 ± 0.05. If n q = 1 is excluded due to its large 
finite-size effect^, the fitting yields z — 1.45 ± 0.07 (solid 
line in Fig^). 

A better approach to determine z is to use the dy- 
namic finite-size scaling theory outlined above. Using 
this method, we obtained z for different values of n q , with 
no resolution function. Such estimates are denoted as z q 
and are shown in FigEl I n our previous worki, we esti- 
mated z as the average value obtained using n q = 1 and 
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FIG. 1: Spin-wave dispersion curve for the average dynamic 
structure factor. 
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FIG. 2: Estimate of the dynamic critical exponent for differ- 
ent n q . Analysis done with L = 30, 36, 48, 60, and 72. 



n q — 2, with a maximum lattice size of L = 60. We now 
have larger L and better statistics, and hence the present 
n q = 1 and n q — 2 correspond to g-values that are closer 
to the asymptotic critical region. Nevertheless, Fig|3 
shows that these values of q are not yet in the asymptotic 
critical region, and a better estimate of z is given by zo, 
obtained by extrapolating z q to the limit q — > 0. We fit- 
ted z q with the function z q — zo + an q + bn q , where zn ; a, 
and b are fitting parameters. Using n q = 1, 2, 7 in the 
fitting (dashed line in Fig0) we find z = 1.48 ±0.02 and 
excluding n q = 1, due to its large finite-size dependence, 
(solid line in FigEJ we find z = 1.50 ± 0.02. Both es- 
timates agree with the theoretical prediction of z = 1.5. 
Dynamic finite-size scaling theory with a small resolu- 
tion function was used to determine z q iteratively. The 
results thus obtained are within a one-cr error bar of the 
respective 5^ = estimates. 
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The longitudinal and transverse components of the dy- 
namic structure factor with respect to the staggered mag- 
netization are shown in Figs^ and 03, respectively, for 
L = 36, T = 0.5T C and q in the [100] direction with 
n q = 3. While the transverse component has a pro- 



0.25 



0.20 



3 



0.15 



0.10 



0.05 



0.00 



0.25 



0.20 - 



3 



0.15 



0.10 



0.05 



0.00 




co/J 



FIG. 3: (a) Longitudinal and (b) transverse components of 
S(q,u>), with respect to the staggered magnetization, with 
5ui = 0.005 at T = 0.5T C . Beside a central peak, the longi- 
tudinal component has two-spin-wave subtraction (indicated 
by solid lines) and addition peaks (dashed lines), correspond- 
ing to 18qi/7r equal to [a] (2,2,0), [b] (1,1,1), [c] (1,1,0), [d] 
(1,0,0), [e] (3,2,0), [f] (3,1,1), [g] (3,1,0), [h] (4,1,1), [i] (2,1,1), 
[k] (1,0,0), [m] (1,1,0), [n] (1,1,1,) and (3,1,0), [o] (3,1,1), [p] 
(2,2,0), [q] (2,1,1) and (3,2,0), [r] (3,2,1) and (4,1,1). 

nounced single spin- wave excitation at to /J « 1.49 and 
intensity ~ 30, as shown in the inset in Fig[3|D, the 
structures on S L (q, oj) have much smaller amplitudes 
and comprise a central peak, and a series of two-spin- 
wave subtraction (a-i) and addition (k-r) peaks. De- 
noting the momentum and frequency of two single spin 
waves as (qi,Wi) and ((^2,^2), the excitations result- 
ing from their addition and subtraction have frequencies 
uj + = loi + u>2 and w_ = \uji — u>2\, respectively, and mo- 
mentum q = qi + q2- For odd values of n q there are no 



two spin-wave peaks at u> = so the central peak seen in 
5 L (q, ui) is presumably due to spin diffusion. 

At T c , S L (q, u) (see Fig0Ji) seems to be predominantly 
diffusive with a central peak that is much more intense 
than the spin- wave peak in 5* T (q, u>) (FigQJ)). These data 
provide clear evidence that the central peak in the aver- 
age S(c[,u>) at T c , seen in experiment and previous sim- 
ulations, but not present in renormalization-group and 
mode-coupling theories, appears only in the longitudinal 
component. Two spin- wave excitations at T c cannot be 
resolved. 



3 

c? 



3 

c? 




FIG. 4: (a) Longitudinal and (b) transverse components of 
S(q,uj), with respect to the staggered magnetization at T c . 

Figure El shows log-log plots of x L (<l) an d X T (<z), the 
integrated intensities of S (q,u) and S T (q,u>), respec- 
tively, as a function of momentum in the [100] direction. 
The momentum dependence of the integrated intensity 
has the form <T X . At T = 0.5T C (FigEH), a linear fit- 
ting in the log-log plane of x L (l) an d X T (q) versus n q for 
L = 36 using n q = 1 to 6 gives \ L ~ l/g°- 80±016 and 
X T ~ l/g 1 - 94±0 06 j whereas if the n q = 1 point is dropped 
the linear fitting yields (solid lines) x L ~ l/g°- 88±0 - 35 
and x T ~ 1/ q lm±o w . These results are in agreement 
with both renormalization-group theory prediction^ and 
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experimental resultsi. We have also tried to fit x L (<z) 
with the Ornstein-Zernike mean- field formula— X L (q) — 
X L (0) K 2 /(q 2 + k 2 ), using x L (0) and K as fitting param- 
eters; however, this expression did not yield a good fit- 
ting. At T c (FigJEJa), our data for L — 60 indicate that 




FIG. 5: Log-log plot of the longitudinal (A) and transverse 
(o) integrated intensity as a function of n q , at (a) T = 0.5T C 
and (b) T = T C . The solid lines are linear fittings using n q = 2 
to n q = 6 for (a), and to n q = 10 for (b). 



X L ~ l/gi-8i±°-°3 and X T ~ l/q l - 92±0M , where n q = 2 to 
10 have been used in the fitting. Finite-size effects on the 
\ow-q divergence exponents of \ L and X T at T c are shown 
in FigEl A linear fitting of these exponents as a function 
of 1/L yields X L ~ l/g 1 - 88 * - 05 and X T ~ l/g 1 - 95 * - 07 
for the thermodynamic limit where L — oo. The dy- 
namic scaling prediction for the static susceptibility at 
T c is§i±£ x = l/? 2 "^ where for the purpose of this com- 
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0.04 ± 0.01. We see that while 



T is consistent with the dynamic scaling prediction, our 



parison we can us> 
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large error bars do not exclude the mean-field behavior of 
X ~ l/? 2 - Our estimate for the divergence of x L at small 
q is slightly less rapid than predicted, but still consistent 
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FIG. 6: Low-g divergence exponents of \ L (A) and x T (°) a t 
T c , as a function of the inverse lattice linear size. The solid 
lines are linear fittings using data for L — 24, 36, 48, 60. 



with it within a two-er error bar. 

Figure [7| shows a log-log plot of the longitudinal and 
transverse components of the characteristic frequency as 
a function of L, for n q — 2 and 6 U = at T c . These 
components are denoted as and w^, respectively, and 
according to dynamic finite-size scaling we have aj^ t = 
L~ zL n L (qL) and w£ = L- zT n T (qL). For n q = 2 (see 
Fig© we obtain z L = 1.48±0.14 and z T = -0.03±0.25, 
where data for L = 36,48, and 60 have been included 
in the analysis. Our data show that the dynamic critical 
exponent for the transverse component is consistent with 
z = 0, indicating that this component is not critical. 
In contrast, the longitudinal component is critical, with 
z w 1.5. 




FIG. 7: Log-log plot of the longitudinal (A) and transverse 
(o) components of the characteristic frequency as a function 
of L with Su, = 0. 
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Separating S'(q, lo) into longitudinal and transverse 
components with respect to the uniform magnetization 
we see that the longitudinal component has a central 
peak, a pronounced spin-wave peak and less intense two- 
spin- wave peaks. In contrast, each such peak in the trans- 
verse component is split into two peaks, with frequencies 
lo l ± Alo, where lo l is the frequency of the corresponding 
peak in the longitudinal component of iS(q, lo). The shift 
Aw in the peak frequencies corresponds to the frequency 
of oscillation of the staggered magnetization. 

IV. CONCLUSIONS 

We have used Monte Carlo and spin dynamics simu- 
lations to study the dynamic behavior of the isotropic 
Heisenberg antiferromagnet on the simple cubic lattice. 
When we use the same range of momentum q as probed 
by experiment, the dynamic critical exponent obtained 7 
is in good agreement with its experimental value, which is 
slightly lower than theoretical predictions. In our present 
work we have used a larger lattice size, and thus smaller 
values of q, in addition to obtaining better statistics and 
extrapolating finite q results to the limit q — 0. This al- 
lowed us to study systematic changes as we approach 
the asymptotic critical region and our improved esti- 
mate (i.e. with systematic errors largely eliminated) is 
z = 1.49 ±0.03, in good agreement with the renormaliza- 
tion group theory and dynamic scaling predictions. Pre- 
sumably the values of q used in the experiment^ and in 
our previous simulations were not in the true asymptotic 
critical region, resulting in a slightly lower estimate of 
z. This should serve as a warning for future simulational 
and experimental probes of dynamic critical behavior. 

Longitudinal and transverse components of S(q, lo) 
with respect to the staggered magnetization are inves- 
tigated separately. This required rotation of the frame 
of reference after each integration step because the stag- 



gered magnetization is not a conserved quantity. Be- 
low T c , the transverse component S T (q,co) has a pro- 
nounced spin-wave peak, whereas the longitudinal com- 
ponent S L (q, lo) is dominated by two-spin- wave addition 
and subtraction peaks, and a central peak presumably 
due to spin diffusion. These results are consistent with 
theory 2 and experiment^, both of which have shown that 
the transverse spin fluctuations are propagating, dom- 
inated by spin-waves, whereas the quasielastic peak is 
due to longitudinal fluctuations. At T c , 5 L (q, lo) has a 
central peak that is much more intense than the spin- 
wave peak in S T (q, lo) and no central peak was seen in 
5 T (q, lo). Explaining the appearance of a central peak in 
5 L (q, lo) at T c remains a challenge for theory. We have 
also seen that while 5 T (q, lo) is not critical, S L (q,ui) is 
critical and it has a dynamic critical exponent z w 1.5. 

These findings further support our earlier conclusion 
that a simple, nearest-neighbor, isotropic Heisenberg 
model describes the behavior of RbMnF3 quite well. The 
only limitations in the agreement appear to be at T c , and 
even there all qualitative features and dynamic exponent 
are faithfully reproduced. 

Below T c our results for the longitudinal and trans- 
verse components of the integrated intensities of 5(q, lo) 
are consistent with renormalization-group theory predic- 
tions, and not with the mean-field one. In contrast, at 
T c , while the integrated intensities are consistent with the 
renormalization-group theory prediction, the large error 
bars do not allow us to exclude the mean-field behavior. 
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